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Abstract 

In  this  paper,  we  construct  a  framework  for  modeling  hysteresis  and  constitutive  nonlinearities  in 
ferroelectric  compounds  based  on  energy  analysis  at  mesoscopic  scales  in  combination  with  stochastic 
homogenization  techniques  to  construct  macroscopic  models.  In  the  first  step  of  the  development, 
previous  analysis  is  used  to  construct  Helmholtz  and  Gibbs  energy  relations  at  the  lattice  level.  This 
provides  local  polarization  relations  which  can  be  extrapolated  to  provide  constitutive  models  for 
certain  homogeneous,  single  crystal  compounds.  To  incorporate  material  and  field  nonhomogeneities, 
as  well  as  the  effects  of  polycrystallinity,  certain  parameters  in  the  local  models  are  assumed  to  be 
manifestations  of  underlying  distributions  having  densities  which  must  be  identified  for  a  given 
compound.  Two  techniques  for  estimating  the  unknown  densities  are  presented,  and  the  accuracy 
of  the  resulting  model  is  illustrated  for  both  symmetric  major  loops  and  biased  minor  loops  through 
fits  and  predictions  with  experimental  PZT4  and  PZT5H  data. 
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1  Introduction 


The  use  of  ferroelectric  compounds  as  actuators  and  sensors  in  high  performance  control  applications 
has  burgeoned  in  recent  years  due  to  the  high  set  point  accuracy  and  broadband  capabilities  of  the 
materials.  For  example,  piezoceramic  rods  or  shells  are  employed  as  positioning  mechanisms  in 
all  present  atomic  force  microscope  (AFM)  and  scanning  tunneling  microscope  (STM)  designs  to 
achieve  the  nanoscale  tolerances  required  by  the  devices  [8,  26,  27]  whereas  relaxor  ferroelectrics 
have  received  significant  attention  as  potential  sonar  transducers  [7,  11]  due  to  the  large  power 
density  to  weight  ratios  exhibited  by  the  materials.  However,  these  advantages  are  accompanied 
by  several  nonlinear  material  attributes  —  hysteresis  is  present  in  the  piezoceramic  materials  at 
all  drive  levels  as  shown  in  Figure  1  whereas  relaxor  ferroelectric  compounds  exhibit  constitutive 
nonlinearities  and  significant  temperature-dependency  for  high  drive  level  operation  in  anhysteretic 
regimes.  These  nonlinear  mechanisms  must  be  accommodated  in  some  manner  to  achieve  the  novel 
control  capabilities  provided  by  these  materials  in  high  performance  applications. 

In  many  cases,  feedback  control  loops  can  be  designed  to  minimize  nonlinear  effects  and  this  has 
led  to  the  successful  use  of  piezoceramic  and  relaxor  ferroelectric  transducers  for  a  broad  range  of 
applications.  However,  this  technique  can  be  difficult  to  implement  at  moderate  to  high  frequencies 
where  increasing  noise-to-signal  ratios  and  diminishing  characteristics  of  high  pass  filters  reduce  the 
effectiveness  of  control  designs  based  solely  on  linear  feedback  laws.  The  use  of  charge  or  current 
controlled  amplifiers  can  also  mitigate  the  majority  of  hysteresis  and  constitutive  nonlinearities  ex¬ 
hibited  by  ferroelectric  materials  [17,  18].  This  can  be  expensive,  however,  when  compared  with 
the  more  common  voltage-controlled  amplifiers,  and  current-controlled  amplifiers  cannot  maintain 
DC  offsets  as  required  by  numerous  applications  —  e.g.,  the  x-state  in  an  AFM  must  hold  a  fixed 
position  during  sweeps  in  the  y-stage.  This  motivates  the  development  of  model-based  frameworks 
which  are  sufficiently  accurate  to  quantify  field,  temperature,  stress,  and  rate-dependent  attributes 
of  ferroelectric  materials  and  sufficiently  efficient  to  permit  real-time  implementation. 

The  PZT5H  data,  plotted  in  Figure  1,  illustrates  a  number  of  the  attributes  which  comprehensive 
models  must  encompass  for  low  frequency,  fixed-temperature,  fixed-stress  regimes.  At  low  field  levels, 
the  materials  exhibit  a  quadratic  Raleigh  loop  behavior  whereas  they  exhibit  saturation  nonlinearities 
and  significant  hysteresis  at  large  field  levels.  In  quasistatic  and  low  frequency  regimes,  models  must 


Figure  1:  Quasistatic  PZT5H  data  collected  at  0.2  Hz  including  a  symmetric  major  loop,  a  symmetric 
Rayleigh  loop,  and  biased  minor  loops. 
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guarantee  closure  of  biased  minor  loops  and  enforce  the  ‘deletion’  or  ‘wiping  out’  property  as  minor 
loops  are  exited  while  including  mechanisms  to  address  accommodation  and  after-effects  (relation 
phenomena).  Additionally,  the  materials  exhibit  stress,  frequency  and  temperature-dependencies 
which  we  address  only  peripherally  in  subsequent  discussion. 

As  detailed  in  [30,  31,  32],  there  exist  a  broad  range  of  techniques  for  modeling  hysteresis  in 
ferroelectric  compounds  including  microscopic  models  [23]  and  both  energy-based  and  purely  phe¬ 
nomenological  macroscopic  models.  We  focus  on  the  latter  category  to  accommodate  the  emphasis 
on  transducer  and  control  design. 

Representative  examples  of  energy-based  modeling  frameworks  are  the  domain  wall  theory  pre¬ 
sented  in  [30,  31],  the  theory  of  Chen  and  Lynch  [4]  and  Huang  and  Tiersten  [12],  and  the  homogenized 
free  energy  theory  [32] .  All  of  these  approaches  have  both  strengths  and  weaknesses  depending  upon 
the  operating  regimes.  For  example,  the  domain  wall  models  are  highly  efficient  for  quantifying 
symmetric  major  loops  but  require  significant  extension  to  guarantee  minor  loop  closure  without  a 
priori  knowledge  of  turning  points  which,  unfortunately,  will  be  the  case  in  feedback  control  appli¬ 
cations.  The  purely  phenomenological  models  include  Preisach  models  which  were  developed  in  the 
context  of  magnetic  materials  [24]  and  have  been  widely  employed  for  a  vast  range  of  applications 
[2,  3]  including  ferroelectric  materials  [10,  25].  The  advantage  of  the  Preisach  approach  lies  in  its 
generality  and  rigorous  mathematical  foundation.  It  has  the  disadvantage  that  this  generality  often 
precludes  the  use  of  physical  measurements  to  identify  parameters  or  update  parameters  to  accom¬ 
modate  changing  operating  conditions.  Furthermore,  the  classical  Preisach  theory  must  be  modified 
in  the  manner  described  in  [9]  for  magnetic  materials  to  accommodate  measured  reversible  effects, 
noncongruency,  and  temperature,  stress  and  rate-dependencies. 

The  theory  presented  here  combines  aspects  of  the  homogenized  free  energy  theory  of  [32]  and 
Preisach  models  posed  in  terms  of  general  densities  or  measures.  In  the  theory  of  [32],  Helmholtz 
and  Gibbs  energy  relations  are  constructed  at  the  lattice,  or  mesoscopic,  level  and  Boltzmann  theory 
is  employed  to  construct  local  polarization  relations.  To  construct  macroscopic  models,  it  is  assumed 
that  local  coercive  and  effective  fields  are  manifestations  of  underlying  lognormal  and  normal  den¬ 
sities.  Stochastic  homogenization  in  this  manner  provides  constitutive  models  for  the  E-P  relation 
which  quantify  hysteresis  and  certain  material  nonlinearities  exhibited  by  ferroelectric  materials. 
However,  the  models  have  limited  accuracy  for  certain  materials  due  to  the  restrictions  imposed 
through  the  a  priori  assumption  of  lognormal  and  normal  densities. 

In  this  paper,  we  reformulate  the  modeling  framework  in  terms  of  general  densities  to  be  identified 
through  a  fit  to  measured  data  from  a  given  compound.  We  also  present  three  identification  strategies 
and  illustrate  the  performance  of  the  model  through  comparison  and  prediction  of  PZT4  and  PZT5H 
data. 

The  formulation  of  Preisach  models  in  terms  of  stochastically  homogenized  free  energy  mod¬ 
els  was  first  proposed  in  [33]  and  the  proposed  framework  extends  this  analysis.  The  construction 
of  hysterons  or  kernels  from  energy  principles  has  the  advantage  that  it  includes  reversible  effects 
directly  and  incorporates  certain  temperature  and  rate  effects  in  the  basis  rather  than  in  the  densi¬ 
ties  or  parameters  as  is  the  case  for  classical  Preisach  models.  Stochastic  homogenization  through 
formulation  in  terms  of  general  densities,  as  motivated  by  Preisach  models,  provides  models  with  sig¬ 
nificant  accuracy  and  ffexibility.  The  resulting  models  guarantee  the  closure  of  biased  minor  loops  in 
quasistatic  or  low  frequency  conditions  but  include  mechanisms  which  characterize  accommodation 
(reptation)  and  after-effect  (relaxation)  phenomena.  Hence  the  framework  incorporates  a  number 
of  the  phenomena  exhibited  by  ferroelectric  materials  while  providing  the  efficiency  required  for 
eventual  implementation. 

From  the  perspective  of  model-based  control  design,  construction  in  terms  of  general  densities 
provides  the  important  advantage  that  models  exhibit  a  linear  dependence  on  the  parameters  as 
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compared  with  the  nonlinear  dependence  which  results  from  the  a  priori  choice  of  lognormal  and 
normal  densities.  This  permits  consideration  of  a  wide  range  of  linear  adaptive  identification  and 
control  strategies  in  the  former  case  whereas  nonlinear  adaptive  identification  strategies  of  the  type 
developed  in  [22]  are  required  in  the  latter  case  and  few  options  exist  for  adaptive  control. 

2  Free  Energy  Model  for  Polycrystalline  Compounds 

We  summarize  here  the  hysteresis  models  developed  in  [32]  for  nonhonrogeneous,  polycrystalline 
ferroelectric  compounds.  In  the  first  step  of  the  model  development,  Helmholtz  and  Gibbs  energy 
relations  are  constructed  at  the  lattice  level  to  quantify  the  balance  between  the  internal  and  elec¬ 
trostatic  energies.  In  the  absence  of  thermally  induced  relaxation  mechanisms,  minima  of  the  Gibbs 
relations  provides  a  mesoscopic  polarization  model  whereas  a  balance  of  the  thermal  and  Gibbs  energy 
through  Boltzmann  principles  is  employed  for  regimes  in  which  relaxation  or  rate-dependent  phe¬ 
nomena  are  significant.  For  homogeneous,  single  crystal  compounds,  the  mesoscopic  models  derived 
through  either  set  of  assumptions  can  be  applied  throughout  the  material  to  provide  a  macroscopic 
model  quantifying  the  hysteresis  and  constitutive  nonlinearities  inherent  to  ferroelectric  materials. 
For  nonhonrogeneous,  polycrystalline  materials,  variations  in  the  lattice  are  incorporated  by  assuming 
that  certain  parameters  in  the  mesoscopic  Gibbs  relations  are  manifestations  of  underlying  distri¬ 
butions  rather  than  constant  values.  In  the  models  derived  in  [32]  and  summarized  in  this  section, 
we  make  the  a  priori  assumption  that  the  local  coercive  fields  Ec  are  distributed  with  a  lognormal 
density  whereas  effective  fields  are  assumed  to  be  normally  distributed  about  the  applied  field  E. 
The  first  assumption  enforces  the  positivity  of  Ec  whereas  the  second  is  based  on  the  tenet  that  local 
dipole  interactions  satisfy  the  central  limit  theorem.  As  illustrated  in  the  examples  of  Section  4  and 
[32],  the  resulting  macroscopic  models  accurately  quantify  the  hysteresis  and  constitutive  nonlinear¬ 
ities  for  PZT5A  and  PZT5H  in  symmetric  drive  regimes  but  have  limited  accuracy  for  hard  PZT4 
compounds  or  drive  regimes  yielding  multiple,  biased  minor  loops. 

In  Section  3,  we  relax  the  assumption  that  Ec  and  Ejt  have  lognormal  and  normal  densities  and 
consider  instead  the  formulation  of  the  model  in  terms  of  general  densities  v  to  be  identified  through 
a  fit  to  data  for  a  given  material.  While  this  requires  additional  complexity  when  constructing  the 
model,  it  yields  superior  accuracy  when  characterizing  general  ferroelectric  compounds  in  complex 
drive  regimes. 

2.1  Mesoscopic  Polarization  Models 

For  fixed  temperature  regimes,  it  is  illustrated  in  [32]  that  an  appropriate  construct  for  the  Helmholtz 
energy  at  the  lattice  level  is 


4(P) 


'  \r] (P  +  PR )2  ,  P<-Pi 

<  HP-Pr)2  ,P>Pi 

k  Hpi  -  pn)  (w  -  P*)  ’  lpl  <  pi- 


(1) 


Here  Pj  and  Pr  respectively  denote  the  inflection  point  in  the  C\  energy  relation  and  the  positive 
point  at  which  a  minimum  of  ^  occurs  as  depicted  in  Figure  2.  The  relation  (1)  quantifies  the  internal 
energy  at  the  lattice  level  in  the  absence  of  applied  fields. 

To  incorporate  the  work  due  to  applied  fields  E,  we  also  consider  the  Gibbs  energy  relation 


G(E,P)  =  -iP(P)  -  EP 


(2) 
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¥(P)=G(0,P) 


Figure  2:  Helmholtz  energy  ^  and  Gibbs  energy  G  for  increasing  field  E  (E2  >  E\  >0).  (b)  Depen¬ 
dence  of  the  local  average  magnetization  P  on  the  field  in  the  absence  of  thermal  activation. 


where  the  latter  term  incorporates  the  potential  energy  U  =  — p  •  E  of  a  dipole  p  in  a  uniaxial  field. 
In  the  absence  of  an  applied  stress  <7,  the  Gibbs  relation  (2)  characterizes  the  energy  landscape  in 
homogeneous,  polycrystalline  materials. 

To  quantify  the  local  average  polarization  P  depicted  in  Figure  2(b),  we  consider  two  cases: 
(i)  drive  regimes  in  which  relaxation  mechanisms  are  negligible,  and  (ii)  regimes  in  which  thermally- 
induced  relaxation  is  sufficiently  prominent  to  warrant  inclusion  in  the  model. 

In  the  absence  of  thermal  relaxation,  the  necessary  conditions  fp  =  0,  >  0  are  invoked  to 

yield  the  general  relation 

P  =  -E  +  5Pr  (3) 

V 

for  the  local  average  polarization.  Here  5  =  1  for  positively  oriented  dipoles  and  5  =  —  1  for 
negative  orientations.  The  expression  (3)  illustrates  that  Pr  can  be  interpreted  as  a  local  or  lattice- 
level  remanence  value  whereas  i]  is  the  reciprocal  slope  after  switching.  To  quantify  the  dipole 
orientations  formally  indicated  by  <5,  we  employ  the  Preisach  notation  (e.g.,  see  [2]) 


[P(E;ECim)=  { 


[P(E-Ec,Om  ,  r(t)  = 


1,-Pr 
l  f  +  Pr 


,  r(t )  /  0  and  F(maxr(i))  =  —  Ec 
,  r(t)  /  0  and  P(maxr(i))  =  Ec. 


Here 


r  E 


[p(E-Ec,om  =  i 


f  -  Pr  ,  E( 0)  <  -Ec 


,  —  Ec  <  E(0)  <  Ec 


If  +  Pr  ,  E(  0)  >  Ec 
denotes  the  initial  dipole  distribution  and  transition  times  are  designated  by 

r(t)  =  {t  <E  (0,  Tf }  |  E(t.)  =  -Ec  or  E(t )  =  Ec} 


(4) 

(5) 

(6) 


where  Tf  denotes  the  final  time  under  consideration.  The  dependence  of  the  kernel  P  on  the  local 
coercive  field  Ec  is  indicated  as  a  prelude  to  the  discussion  in  Section  2.2  where  Ec  is  assumed  to  be 
distributed. 
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For  operating  regimes  in  which  relaxation  mechanisms  are  significant,  it  is  necessary  to  balance 
the  Gibbs  energy  G  with  the  relative  thermal  energy  kT /V  —  where  k  denotes  Boltzmann’s  constant 
—  through  the  Boltzmann  relation 

H(G)  =  Ce~GV/kT.  (7) 

It  is  observed  that  the  probability  p  of  obtaining  an  energy  level  G  is  increased  when  the  magnitude 
of  kT /V  approaches  that  of  G.  The  constant  C  is  chosen  to  ensure  a  probability  of  one  for  integration 
over  all  possible  dipole  orientations. 

If  we  let  x+  and  X-  denote  the  fractions  of  dipoles  having  positive  and  negative  orientations,  and 
let  (P+)  and  (P_)  denote  the  average  expected  polarizations  associated  with  the  two  orientations, 
the  local  average  polarization  at  the  lattice  level  is 


P  =  x+  ( P+ }  +  X-  ( P- ) . 


(8) 


Since  the  expected  polarization  values  are  obtained  by  integrating  the  product  Pp(G(P))  over  all 
admissible  configurations,  it  follows  that 


(P+) 


pe-G(E,P,T)V/kTdp 


/  e-G(E,P,T)V/kTdp 

I  Pi 


r~Pi 


pe-G(E,P)V/kTdp 


(P-)  = 


r~Pi 


-G(E,P)V/kTdp 


(9) 


The  denominator  results  from  the  evaluation  of  the  integration  constant  C  whereas  it  is  illustrated  in 
[32]  that  the  use  of  the  inflection  points  ± Pi,  to  simplify  evaluation  of  the  integrals,  can  be  justified 
though  either  asymptotic  analysis  or  energy  arguments. 

Debye  arguments  yield  the  differential  equations 


X+  =  —p-i I _ x_|_  +p _ |_x_ 

X-  =  —p _ yX-  +  P- 1 _ X+ 


(10) 


quantifying  the  evolution  of  the  respective  dipole  fractions.  For  implementation,  the  relations  (10) 
can  be  simplified  to  the  single  differential  equation 


x+  =  -P+-X+  +P-+(  1  -  x+) 


(11) 


through  the  identity  x+  +  x-  =  1. 

The  likelihoods  of  switching  from  positive  to  negative  orientations,  or  vice  versa,  are  respectively 
quantified  by 


rPi+e 


e-G(E,PI,T)V/KTdp 


P+-  = 


I Pi 


T(T) 


e-G(E,P,T)V/kT dp 


'Pi 


P-  + 


1 

nf) 


G{E-PI,T)V/kTdp 


G(E,P,T)V/kT  dp 


(12) 


where  e  is  taken  to  be  a  small  positive  constant.  The  relaxation  term  T  quantifies  the  frequency 
at  which  jumps  are  attempted  whereas  the  remainder  of  the  definition  characterizes  the  probability 
of  achieving  the  energy  required  to  exit  respective  potential  wells.  It  is  detailed  in  [32]  that  this 
probability  increases  when  the  relative  thermal  energy  kT /V  approaches  the  Gibbs  energy  G.  For 
regimes  in  which  kT /V  «  G,  it  is  also  illustrated  that  the  average  local  polarization  relation  (8) 
asymptotically  limits  to  the  expressions  (3)  or  (4)  employed  for  negligible  thermal  relaxation. 
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The  polarization  relations  (4)  and  (8)  are  derived  at  the  lattice  level  and  hence  are  mesoscopic 
in  nature.  For  single  crystal,  homogeneous  materials  with  uniform  effective  fields,  the  relations  can 
be  extrapolated  throughout  the  material  to  yield  constitutive  models  quantifying  the  hysteretic  and 
nonlinear  dependence  of  the  bulk  polarization  P  on  input  fields  E.  Macroscopic  models  constructed 
in  this  manner  retain  the  sharp  transition  through  the  coercive  point  exhibited  by  the  kernels  (4)  or 
(8).  While  this  models  the  physical  behavior  of  certain  single  crystal  compounds  such  as  BaTiC>3  — 
e.g.,  see  page  76  of  [21]  —  the  transition  is  steeper  than  that  typically  exhibited  by  nonhomogeneous, 
polycrystalline  compounds.  To  provide  macroscopic  E-P  models  appropriate  for  this  broader  class 
of  compounds,  stochastic  homogenization  techniques  utilizing  the  mesoscopic  relations  (4)  or  (8)  as 
kernels  are  considered  in  the  next  section. 

2.2  Macroscopic  Polarization  Models 

The  local  average  polarization  relations  (4)  and  (8)  were  derived  under  the  assumption  of  a  uniform 
lattice  and  constant  effective  field  Ee  =  E.  Hence  they  do  not  incorporate  variability  in  the  lattice  due 
to  nonunifornr  grain  structure,  material  nonhomogeneities  and  defects,  and  variable  stress  fields  nor 
do  they  include  the  effects  of  nonunifornr  interaction  or  effective  fields.  These  effects  can  in  theory  be 
incorporated  through  nricronrechanical  models  utilizing  extended  energy  relations  [5,  6,  16].  However, 
this  produces  models  whose  complexity  precludes  bulk  material  characterization,  transducer  design  or 
model-based  control  design.  Alternatively,  one  can  assume  that  certain  parameters  in  the  mesoscopic 
energy  relations  are  manifestations  of  underlying  stochastic  distributions  rather  than  constant  values 
as  posited  for  single  crystals  having  uniform  lattices.  Stochastic  homogenization  in  this  manner 
produces  macroscopic  models  which  retain  energy  characteristics  but  are  sufficiently  low-order  to 
permit  implementation. 

As  noted  in  [32],  variability  in  the  lattice  can  be  incorporated  by  considering  the  parameters 
Pr,Pi  or 

Ec  =  v{Pr  ~  Pi) 

to  be  distributed,  whereas  variable  interaction  fields  can  be  incorporated  by  considering  effective 
fields  Ee  to  be  distributed  about  the  applied  field  E.  In  this  section,  we  summarize  the  models  of 
[32]  based  on  a  priori  choices  for  the  distributions  whereas  in  Section  3,  we  consider  the  development 
of  models  having  general  densities  estimated  through  fits  to  experimental  data. 

As  illustrated  in  Figure  2,  the  coercive  field  Ec  is  nonnegative  which  constitutes  a  constraint  on 
the  underlying  density.  Two  a  priori  choices  which  enforce  this  constraint  are  normal  distributions 
truncated  to  include  only  nonnegative  values  and  lognormal  distributions.  The  latter  are  employed 
in  [32]  which  yield  the  density 

Vl(Ec)  =  Cie-MEcfEc)/2c]*  (13) 

where  ci,  c  and  Ec  are  positive  constants.  It  is  illustrated  in  [9]  that  if  c  is  small  compared  with  Ec, 
the  mean  and  variance  for  the  distribution  have  the  approximate  values 

(Ec)  «  Ec  ,  a^2Ecc .  (14) 

It  is  illustrated  in  the  examples  of  Section  4  that  (14)  provides  initial  estimates  for  the  parameters 
Ec  and  c  based  on  measured  properties  of  the  data.  We  note  that  lognormal  distributions  of  the  form 
(13)  are  employed  when  quantifying  the  distribution  of  the  critical  field  parameter  Hk  in  Preisach 
models  used  to  model  hard  magnetic  materials  [9]  which  provides  additional  motivation  for  their 
consideration  in  the  present  setting. 

Secondly,  we  consider  variations  in  the  effective  field  at  the  lattice  level.  As  noted  in  [1,  20,  30,  31], 
the  applied  field  E  is  augmented  by  an  interaction  field  E)  due  to  neighboring  dipoles  as  well  as 
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certain  electromechanical  interactions.  While  microelectric  energy  analysis  can  quantify  some  of 
these  mechanisms,  the  required  complexity  of  subsequent  models  precludes  transducer  design  or  real 
time  control  implementation.  Alternatively,  we  make  the  a  priori  assumption  that  effective  fields 
Ee  =  E  +  Ei  are  distributed  about  the  applied  field  with  an  underlying  density  characteristic  of  the 
ferroelectric  material  under  consideration.  The  models  employed  in  [32]  are  based  on  the  assumption 
that  Ee  is  normally  distributed  about  E  with  the  density 

u2(Ee-1E)  =  c2e-(E-E^/2b2.  (15) 

The  magnitude  of  the  variance  b 2  dictates  the  degree  to  which  switching  occurs  prior  to  the  remanence 
polarization.  Large  values  of  b2  produce  models  with  significant  switching  as  the  applied  field  E  is 
reduced  to  zero  with  correspondingly  large  changes  in  the  slope  of  the  hysteresis  curve  whereas  small 
values  of  b2  are  employed  when  modeling  the  hysteretic  response  of  materials  exhibiting  nearly  linear 
E-P  behavior  at  remanence. 

The  resulting  polarization  model  is 

/*oo  />oo  _ 

[P(E)\(t)  =  C  /  {PiE  +  EuE^Ome-^/^e-^^^/^dEidE,  (16) 

J  0  J — oo 

where  the  kernel  P  is  specified  by  (3),  (4)  or  (8)  and  £  denotes  the  initial  distribution  of  dipoles. 


2.3  Discrete  Macroscopic  Polarization  Model 

To  implement  the  model  (16),  it  is  necessary  to  approximate  the  integrals.  As  detailed  in  [32],  this 
can  be  accomplished  either  by  employing  Gaussian  quadrature  routines  constructed  for  infinite  or 
semi-infinite  domains  or  by  exploiting  the  decay  exhibited  by  the  kernels  to  truncate  the  domains 
and  employ  Gauss-Legendre  quadrature  rules.  Both  strategies  yield  discrete  models  of  the  form 


N,,  Nj 


[P(E)](t)  =  C^2^2[P(E  +  Et 


i',ECi,€i)](t)e  Eii/be  MWWviWi 


i=l  j=  1 


(17) 


where  E^ ,  ECi  denote  the  abscissas  associated  with  respective  quadrature  formulae  and  Wj  are  the 
respective  weights  —  e.g.,  see  pages  698-699  of  [37]. 

For  the  characterization  examples  presented  in  Section  4,  as  well  as  in  [32],  we  employ  a  composite, 
4  point  Gauss-Legendre  quadrature  rule  constructed  for  truncated  domains,  and  we  illustrate  that 
approach  here.  Consider  first  the  approximation  of  the  effective  field  integral.  For  a  given  threshold 
e  and  index  Nq,  we  let  [—L,L\  denote  the  interval  where  \e~Ei  \  >  e  and  consider  the  partition 
hq  =  —L  +  qh,h  =  2 L/Nq.  On  each  subinterval  [hq-i,hq],  the  quadrature  points  and  weights  are 
specified  to  be 


E„ 


1 ql 


E, 


1q2 


E, 


Iqi 


E, 


hq—  1  +  h 
hq—  1  +  h 
hq—  1  +  h 
hq—  1  +  h 


\J  15+2\/30 

„  .  _  4 9h 

2\/35 

ql  12(18+^30) 

\J  15— 2\/30 

_  49/i 

2\/35 

’  12(18-v/30) 

V  15-2%/30 

„..  _  49/i 

2\/35 

93  12(18— \/30) 

Vl5+2V30 

.  —  49/i 

2^35 

’  Wq4  ~  12(18+^30) 

(18) 
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Similar  analysis  yields  the  quadrature  points  ECpk  and  weights  vpk  for  the  approximation  of  the 
coercive  field  integral  on  a  grid  with  Np  subintervals.  In  this  case,  TVj  =  ANp,  Nj  =  4:Nq,  and  the 
formulation  (17)  has  the  specific  form 


[P(E)](t)  =  C 


Np  4 

EE 

p=i  k=i 


Nq  4 

E  E^V  +  ^  ECpk,^)](t)e^/2b2 

q=  1  1=1 


Wqt 


-[\n(ECpk/Ec)/2c}< 


Vpk  ■  (19) 


Details  regarding  the  discrete  formulation  model  and  a  highly  efficient  implementation  algorithm  are 
provided  in  [32], 

We  note  that  for  the  linear  kernel  P  given  by  (3)  or  (4),  construction  of  the  model  using  this 
a  priori  choice  of  normal  and  lognormal  distributions  entails  the  identification  of  6  parameters.  The 
parameters  Pr  and  rj  arise  in  the  definition  of  the  kernel  whereas  Ec,c,b  and  C  are  associated  with 
the  construction  of  the  density  functions  u\  and  v2.  We  note  that  for  this  choice  of  kernel,  Pr  and 
C  both  scale  final  polarization  values  and  hence  can  be  combined  into  a  single  effective  parameter. 
Thus  for  materials  in  which  the  normal  and  lognormal  densities  provide  sufficient  accuracy,  the  low 
number  and  physical  nature  of  parameters  makes  the  model  highly  efficient  to  construct  and  update. 


3  Polarization  Model  with  General  Densities 

Whereas  certain  physical  arguments  motivate  consideration  of  the  normal  and  lognormal  distributions 
for  the  effective  and  coercive  fields,  we  rely  primarily  on  their  mathematical  attributes  including  the 
easily  quantified  relations  between  properties  of  the  data  —  including  the  coercivity  and  degree  of 
pre-remanence  switching  —  and  mean  and  variance  properties  of  the  densities.  Moreover,  as  long  as 
the  densities  enforce  positive  arguments  for  the  coercive  field  and  exhibit  certain  decay  properties, 
there  is  no  physical  reason  why  normal  and  lognormal  functions  are  inherently  preferable  to  general 
densities  which  provide  substantially  greater  flexibility  for  model  construction.  We  thus  consider  the 
formulation  of  the  polarization  model  in  terms  of  general  densities  and  indicate  techniques  suitable  for 
estimating  the  densities.  The  attributes  of  models  constructed  using  general  densities  are  compared 
with  those  based  on  the  a  priori  assumption  of  normal  and  lognormal  distributions  in  Section  4. 


3.1  Polarization  Model 


We  consider  general  densities  u\  and  v2  which  satisfy  the  conditions 

(i)  v\  (x)  defined  for  x  >  0, 

(ii)  v2(-x)  =  v2(x),  (20) 

(hi)  Mx)|  <Cle~aix  ,  \is2(x)\  <  c2e~a2X 


for  nonnegative  ci,ai,C2,a2.  The  restricted  domain  in  (i)  reflects  the  fact  that  the  coercive  field 
Ec  is  positive  whereas  the  symmetry  enforced  in  the  effective  field  through  (ii)  yields  the  symmetry 
observed  in  low- field  Rayleigh  loops.  Hypothesis  (iii)  incorporates  the  physical  observation  that 
the  coercive  and  interactions  fields  decay  as  a  function  of  distance  and  guarantees  that  integration 
against  the  piecewise  linear  kernel  yields  finite  polarization  values. 

The  general  polarization  model  can  then  be  formulated  as 


[p(E)m  = 


/*oo  /*oo 

'  0  J — oo 

roo  /»oo 

1 0  J — oo 


vi{Ec)v2(Ei){P(E  +  Ep  Ec,  £)](f)  dEi  dEc 
e{Ec,  Ei)  [P(E  +  Ep  Ec,£)](t)  dE{  dEc 


(21) 


where  P  is  specified  in  (3),  (4)  or  (8).  Whereas  formulation  of  the  model  in  terms  of  the  two- 
dimensional  density  v  is  more  general,  retention  of  the  components  u\  and  v2  can  facilitate  subsequent 
implement  ation . 

The  approximation  of  the  integrals  yields  the  discrete  model 

Ni  Nj 

[P(E)](t)  +  (22) 

i=l  j= 1 

where  Vi,Wj  denote  the  quadrature  points  and  ECi,Eij  are  the  abscissas  for  the  given  quadrature 
rule. 

To  illustrate  attributes  of  the  discrete  model  (22),  we  compare  it  with  the  model  (17)  derived 
under  the  a  priori  assumption  that 


v'i(Ec)  =  Cie-K£c/sc)/2c]’ 
v2  =  c2e~E^b. 


(23) 


The  construction  of  (17)  requires  the  identification  of  the  5  parameters  rj,C,Ec,c  and  b  once  Pr  is 
combined  with  C.  The  construction  of  the  model  (22)  in  terms  of  v\  and  v2  requires  the  identification 
of  the  Ni  +  Nj  +  1  parameters  [v\ (EC1),  •  •  • ,  (ECn,  )],  \v2(Eif),  ■  ■  ■ ,  v2{EjN  )]  and  ip  If  one  employs 
the  product  density  u,  discretization  yields  Nt  ■  Nj  +  1  parameters  since  one  must  identify  all  of  the 
values  v(ECi ,  Eij )  in  this  formulation.  Since  one  typically  requires  Ni  and  Nj  on  the  order  of  20  to  80 
to  achieve  convergence,  the  identification  of  the  general  densities  using  either  the  isolated  densities 
v\  and  v2 ,  or  the  product  density  u,  requires  highly  efficient  optimization  and  parameter  estimation 
techniques. 

The  advantage  gained  through  this  effort  is  the  construction  of  models  which  are  highly  accurate 
for  a  wide  range  of  drive  regimes.  It  is  important  to  note  that  whereas  the  general  density  model 
(23)  may  be  more  expensive  to  construct  than  the  model  (17),  the  implementation  of  the  two  models 
requires  identical  overhead  since  both  simply  employ  multiplication  of  Nt  x  1  and  Nj  x  1  vectors. 
Hence  they  will  be  equally  efficient  for  model-based,  control  design. 


3.2  Estimation  of  Ui  and  v2  through  Constrained  Optimization 

To  estimate  the  N  +  Nj  +  1  parameters  [iq (ECl),- ■  ■  ,ui(ECN  )\,[i'2(Ei1),-  ■  ■  ,u2(EiN  )\  and  ip  we 

_  *  3 

consider  a  least  squares  fit  to  data  (£&,  P*.),  k  =  1,  •  •  ■ ,  Nrj.  The  accuracy  of  resulting  models  will 
be  improved  if  the  data  is  chosen  to  include  all  drive  regimes  under  consideration  and,  in  general, 
inclusion  of  a  highly  varied  set  of  drive  regimes  will  provide  more  comprehensive  characterization 
of  the  densities.  For  example,  identification  of  the  densities  solely  based  on  symmetric  major  loop 
data  will  provide  a  model  which  has  moderate  accuracy  when  predicting  biased  minor  loops  whereas 
identification  using  data  that  includes  biased  drive  level  data  will  provide  a  model  with  improved 
accuracy  in  these  regimes. 

To  formulate  the  constrained  optimization  problem,  we  define  6  =  [rj,  v\(ECl  ),•■•,  v\{ECn), 
i^2 ( Ej , ) ,  •  •  • ,  u2(EiN  )\  E  ] YiNi+Nj+i  ancj  iet  P(Pfc;  0)  denote  the  modeled  parameter-dependent  polar¬ 
ization  values  given  by  (22)  for  the  7Vfi  input  field  values  E^.  The  constrained  optimization  problem 
can  then  be  formulated  as  follows:  find  0  E  lRNi+N i+l  which  minimizes 

f(0)=1-\\P(Ek]0)-Pk\\2  (24) 
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subject  to 


Oi  >  0. 

Here  ||  •  ||  denotes  the  Euclidean  norm  in  ETVd.  For  the  general  formulation  (24),  optimization  can 
be  performed  using  the  Matlab  routine  fmincom.m,  or  other  optimization  routines  which  enforce  the 
positivity  constraint  for  moderate  values  of  Nt  and  Nj  —  e.g.,  Nt  =  Nj  =  80.  We  note  that  solution  of 
(24)  is  analogous  to  the  technique  employed  in  [34]  for  estimating  parameters  in  discretized  Preisach 
models. 


3.3  Estimation  of  v  —  Direct  Solution  Algorithm 

Alternatively,  one  can  consider  the  more  general  formulation 


[p(E)m  = 


R{EC,  Ei)[P(E  +  Ej]  Ec,  £)](f)  dEj  dEc , 


(25) 


where  v  :  ]R2  — >  ]R  is  a  general  density  function  to  be  estimated  for  a  given  compound,  and  discretized 
relation 

Ni  Nj 

[P(E)](t)  =  YJY,^ECi,Eij)[P{E  +  Eij-ECi,m)viWj.  (26) 

i=  1  3= 1 

This  formulation  has  the  advantage  that  it  yields  a  linear  system  in  terms  of  v.  Furthermore,  it 
yields  a  quadratic  programming  problem  when  the  discrete  model  (26)  is  employed  in  a  least  squares 
setting  but  it  comes  at  the  cost  of  estimating  Nt  ■  Nj  unknowns  as  compared  with  the  Nj  +  Nj  +  1 
variables  required  for  the  nonlinear  constrained  optimization  problem  discussed  in  Section  3.2  —  we 
assume  here  that  rj  can  be  estimated  directly  from  the  slope  ^  of  data  after  switching. 

In  the  absence  of  thermal  relaxation,  evaluation  of  (3)  at  the  measured  input  field  values  Ek 
yields 

Nj  n  r 


P(Ek)  =  EE 

j= i *=i  L 


Ek  +  Ej, 


V 


+  PRS(Ek-  ECi ,  Ej. ) 


v(ECi,Ei.)viWj 


(27) 


for  k  =  1,  •  •  • ,  N(j.  For  clarity,  we  note  the  dependence  of  6  on  the  coercive  field  quadrature  points. 
To  formulate  (27)  as  a  linear  system,  we  define  the  Nj  x  Nj  matrices  Ak  and  to  have  components 


Ek  +  Ej. 
V 


+  PRS(Ek-  ECj ,  Ej. 


VjW 


(28) 


[^] ij  ~  ^{Ea  >  Ejj ) . 


For  N  =  Nj  ■  Nj,  we  define  the  N  x  1  vector  9  and  1  x  N  vector  ak  by 

9  =  vec(<f>)  ,  ak  =  [vec(Afc)]T 


(29) 


where  ‘vec’  denotes  the  vector  concatenation  of  the  respective  matrices.  Additionally,  the  N^  x  1 
vectors  V  and  V  are  defined  componentwise  by 


[P]k  =  P{Ek)  ,  [P]k  =  Pk • 

Finally,  the  x  N  matrix  A  is  defined  row-wise  by 


(30) 


[A\k  =  ak. 
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The  polarization  model  (27)  can  then  be  formulated  as  the  linear  system 


A9  =  V  (31) 

and  the  least  squares  problem  used  to  identify  the  density  values  6  given  measurements  {E^,  P^},  k  = 
1,  •  •  • ,  Nd,  is  the  following:  minimize 


subject  to 


m  =  1-\\ao-v\\2 


Because  the  minimum  is  unaffected  by  a  constant  shift,  this  is  equivalent  to  the  problem 

min  f(6)  ,  f(0)  =  \eT At A6  -  VT A9 
6  2 

subject  to  9j  >  0 ,  j  =  1,  •  •  • ,  IV. 


(32) 


(33) 


To  solve  (33),  we  employ  the  same  technique  used  in  [28]  for  identifying  density  functions  in 
Preisach  models.  Consider  the  singular  value  decomposition 

AT  A  =  USVT,  (34) 

where  U  =  V  since  AT A  is  symmetric,  and  the  N  x  N  diagonal  matrix  S  is  comprised  of  the  singular 
values  of  ATA.  Furthermore,  it  is  noted  that  the  column  vectors  Ui  of  U  satisfy  the  orthogonality 
condition 

ufuj  =  Sij  (35) 

where  5^  denotes  the  Kronecker  delta. 

We  consider  the  case  rank(yl)  =  rank (5')  =  q  <  minjlV,  Nj}  which  arises  when  considering  fine 
discretizations  to  fully  resolve  fine-scale  properties  of  the  hysteretic  response.  For  this  case,  we 
eliminate  the  rows  and  columns  of  S  and  U  corresponding  to  zero  singular  values  to  obtain  the  q  X  q 
matrix  S  and  N  x  q  matrix  U.  If  instability  due  to  ill-conditioning  induced  by  small  singular  values 
is  of  concern,  one  can  alternatively  retain  singular  values  greater  than  a  specified  threshold  e  along 
with  corresponding  components  in  U  —  this  is  the  basis  for  the  truncated  SVD  (TSVD)  techniques. 
Because 

AtA  =  USUT  ,  UTU  =  I, 

the  minimization  problem  (33)  can  be  reformulated  as 

min  f(x)  ,  fix)  =  -xT Sx  —  VT AUx 

*  2  (36) 

subject  to  g(x)  =  Ux  >  0 

where  x  =  UT0  and  hence  0  =  Ux.  The  quadratic  programming  problem  (36)  can  be  solved  using 
the  Matlab  routine  quadprog.m  to  obtain  x*,  and  a  corresponding  solution  9*  to  (33)  is  then  given 
by  9*  =  Ux*. 
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3.4  Estimation  of  v  —  Regularized  Solution 

It  is  shown  in  [29]  that  the  polarization  model  (25)  is  an  integral  equation  having  a  compact  integral 
operator.  Hence  the  inverse  problem  of  determining  v  given  data  measurements  is  ill-posed  both 
with  regard  to  the  existence  of  a  unique  solution  and  the  continuous  dependence  of  u  on  the  data. 
It  is  illustrated  in  the  validation  example  of  Section  4.3  that  this  ill-posedness  is  manifested  in  the 
discretized  formulation  (26)  through  a  loss  in  accuracy  as  discretization  limits  are  increased.  This 
motivates  the  consideration  of  regularized  least  squares  formulations  for  determining  v. 

As  detailed  in  [29],  Tikhonov  regularization  is  imposed  by  replacing  the  least  squares  functional 
(32)  by  the  augmented  functional 

/«(«)  =  1m»-p||2+||I«II2  (37) 

which  shifts  the  spectrum  of  AT A  to  ATA  +  al  to  avoid  the  deleterious  effects  of  small  singular 
values.  The  constrained  minimization  problem  in  this  case  is 

min  fa(0)  ,  f(9)  =  ^\\A9-V\\2  +  |||0||2 

subject  to  9j  >  0,  j  =  1, . . . ,  N. 

Techniques  for  choosing  a  to  avoid  oversnroothing  solutions  as  well  as  a  solution  algorithm  for  (38) 
can  be  found  in  Vogel  [36] . 

4  Model  Validation 

To  illustrate  attributes  of  the  model  and  the  flexibility  provided  by  the  formulation  in  terms  of  gen¬ 
eral  densities,  we  consider  the  characterization  of  PZT5H  and  PZT4.  The  PZT5H  example  illustrates 
the  capability  of  the  model  to  accurately  characterize  biased  minor  loops  under  quasistatic  drive  con¬ 
ditions.  An  important  attribute  of  the  model  is  its  capability  for  ensuring  closure  of  the  loops  under 
such  conditions  in  accordance  with  experimental  behavior  of  the  materials.  The  PZT4  example  illus¬ 
trates  the  predictive  capabilities  for  the  model  for  a  hard  PZT  compound.  In  both  cases,  we  employ 
the  polarization  kernel  (3)  or  (4)  since  relaxation  mechanisms  are  negligible.  For  operating  regimes  in 
which  relaxation  must  be  accommodated,  the  kernel  (8)  would  be  employed  in  an  analogous  manner. 

4.1  PZT4  Characterization 

We  consider  first  the  characterization  of  constitutive  nonlinear  and  hysteresis  exhibited  by  PZT4. 
Data  was  collected  from  a  3.81  cm  x  0.635  cm  x  0.381  cm  wafer  at  200  mHz  to  minimize  frequency 
effects.  The  measured  E-P  relations  corresponding  to  peak  input  voltages  ranging  from  600  V  to 
1800  V  are  plotted  in  Figure  3.  Note  that  the  relation  between  E  and  V  can  be  approximated  by 
E  =  V/h  where  h  =  3.81  x  10-4  is  the  thickness  of  the  wafer.  To  demonstrate  properties  of  the 
model  using  both  the  normal/lognormal  and  general  densities  and  differing  identification  strategies, 
we  consider  four  regimes: 

(i)  normal/lognormal  densities  estimated  using  all  four  data  sets; 

(ii)  general  densities  estimated  using  all  four  data  sets; 

(iii)  normal/lognormal  densities  estimated  using  the  1800  V  data  set; 

(iv)  general  densities  estimated  using  the  1800  V  data  set. 

In  all  cases,  convergence  was  achieved  with  Np  =  Nq  =  20  which  yields  the  parameter  limits  ATt  = 
Nj  =  80. 


12 


Lognormal/Normal  Densities  —  Identification  Using  600  V  -  1800  V  Data 

To  provide  a  baseline  for  comparison,  we  take  u\  and  V2  as  the  lognormal  and  normal  functions 
specified  in  (23)  and  employ  four  data  sets  —  600  V,  1000  V,  1200  V  and  1800  V  —  when  estimating 
the  parameters  Pr,t],  Ec,c,b  and  C.  The  resulting  model  response,  is  compared  with  the  data  in 
Figure  3.  The  coercive  and  effective  field  densities  v\  and  v-2  are  plotted  in  Figure  4.  It  is  observed 
that  the  model  is  reasonably  accurate  throughout  the  drive  range  which  is  to  be  expected  since  all 
four  data  sets  were  employed  when  estimating  parameters.  The  accuracy  is  quantified  by  the  residual 


n  = 


i 

Nd 


E 


k=l 


1/2 


| P(Ek)  -  Pk |2 


(39) 


between  the  modeled  response  and  polarization  data  at  the  specified  input  field  values  which,  as 
summarized  in  Table  1,  is  0.0226. 


(i) 

(ii) 

(hi) 

General  Densities 

0.0176 

0.0383 

0.0113 

Lognormal/Normal  Densities 

0.0226 

0.0533 

0.0164 

Table  1:  Residuals  1Z  given  by  (39).  (i)  Identification  and  residuals  over  four  data  sets  (600  V  -  1800  V 
data),  (ii)  Identification  using  1800  V  data  —  residual  for  all  four  data  sets,  (iii)  Identification  using 
1800  V  data  —  residual  for  1800  V  data. 


Figure  3:  PZT4  data  ( - )  and  model  predictions  ( - )  with  lognormal/normal  densities  for  zq 

and  i>2  and  parameters  estimated  through  a  least  squares  fit  to  all  four  data  sets.  Abscissas:  electric 
field  (MV/m),  ordinates:  polarization  (C/m2). 
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Coercive  Field  (MV/m)  Effective  Field  (MV/m) 


Figure  4:  Lognormal  coercive  field  density  v\  and  effective  field  density  zq  given  by  (23)  with  pa¬ 
rameters  estimated  through  a  least  squares  fit  to  600  V  -  1800  V  data. 

General  Densities  —  Identification  Using  600  V  -  1800  V  Data 

We  now  assume  that  v\  and  zq  are  general  functions  satisfying  the  hypotheses  (20)  and  we 
estimate  the  JVj  +  Nj  +  1  parameters  [zq (ECl),  ■  ■  ■ ,  zq (ECn  )\,  •  •  • ,  zq (EiN  )\  and  rj,  through 

t  j 

the  solution  of  the  constrained  least  squares  problem  (24).  The  resulting  model  fits  are  plotted  in 
Figure  5  and  the  estimated  densities  are  plotted  in  Figure  6.  A  comparison  between  Figures  3  and  5 
indicates  that  the  model  with  general  densities  is  more  accurate  as  quantified  by  the  residual  of  0.0176 
which  is  approximately  60%  that  obtained  using  the  a  priori  choice  of  lognormal  and  normal  densities 


Figure  5:  PZT4  data  ( - )  and  model  predictions  ( - )  obtained  with  the  general  densities  zq 

and  z/2  estimated  through  a  constrained  least  squares  fit  to  all  four  data  sets.  Abscissas:  electric  field 
(MV/m),  ordinates:  polarization  (C/m2). 
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Figure  6:  Coercive  field  density  u\  and  effective  field  density  U2  estimated  using  all  four  data  sets. 


—  see  Table  1.  A  comparison  between  the  coercive  and  effective  field  densities  and  those  plotted  in 
Figure  4  illustrates  that  the  general  densities  have  a  qualitative  shape  similar  to  the  lognormal  and 
normal  densities  but  are  significantly  less  regular  and  do  not  exhibit  the  same  nronotonicity. 

Lognormal/Normal  Densities  —  Identification  Using  1800  V  Data 

The  third  case  consists  of  identification  of  the  lognormal  and  normal  densities  using  the  1800  V 
data  and  subsequent  model  predictions  at  the  600  V,  1000  V  and  1200  V  levels.  This  yielded  the 
model  responses  plotted  in  Figure  7.  A  comparison  between  these  results  and  those  plotted  in 
Figure  3  illustrates  a  slightly  improved  fit  at  1800  V  but  a  significantly  worse  fit  at  600  V.  This  is 


Figure  7:  PZT4  data  ( - )  and  model  predictions  ( - )  with  lognormal/normal  densities  for  vi 

and  v 2  and  parameters  estimated  through  a  least  squares  fit  to  the  1800  V  data.  Abscissas:  electric 
field  (MV/m),  ordinates:  polarization  (C/m2). 
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reflected  in  the  residual  of  0.0533  which  is  more  than  a  factor  of  two  worse  than  when  all  four  data 
sets  are  employed  for  parameter  estimation.  This  illustrates  that  while  the  model  has  substantial 
predictive  capabilities,  improved  performance  can  be  obtained  by  estimating  parameters  using  data 
from  a  range  of  drive  regimes. 

General  Densities  —  Identification  Using  1800  V  Data 

To  complete  the  repertoire,  we  consider  the  identification  of  the  general  densities  v\  and  v2 
through  solution  of  the  constrained  minimization  problem  (24)  using  only  the  1800  V  data.  A 
comparison  of  the  resulting  model  fits  in  Figure  8  with  those  in  Figure  5  illustrates  an  improved 
and  highly  accurate  fit  to  the  1800  V  data,  reasonable  predictions  at  1000  V  and  1200  V,  and 
only  moderately  accurate  predictions  at  600  V.  As  noted  in  Table  1,  the  residual  of  0.0383  over  all 
four  data  sets  is  larger  than  those  obtained  with  either  general  densities  or  the  a  priori  choice  of 
lognormal/normal  densities  when  all  four  data  sets  are  employed  for  identification  but  is  significantly 
smaller  than  the  lognormal/normal  model  fit  when  identifying  parameters  using  only  the  1800  V  data. 
This  same  tendency  is  noted  in  the  residuals  over  only  the  1800  V  data. 


Figure  8:  PZT4  data  ( - )  and  model  predictions  ( - )  obtained  with  the  general  densities  01 

and  v2  estimated  through  a  constrained  least  squares  fit  to  the  1800  V  data.  Abscissas:  electric  field 
(MV/rn),  ordinates:  polarization  (C/m2). 


4.2  PZT5H  Characterization  —  Individual  Densities  and  v2 

The  second  compound  we  consider  is  PZT5H  which  exhibits  the  hysteretic  and  nonlinear  behavior 
shown  in  Figure  1.  In  this  case,  we  illustrate  the  capability  of  the  model  to  accurately  quantify 
biased  minor  loop  behavior  in  quasistatic  drive  regimes  —  data  was  collected  at  0.2  Hz.  We  again 
consider  four  cases  to  illustrate  the  performance  of  the  model  constructed  with  general  densities  as 
compared  with  the  model  employing  lognormal/normal  densities: 
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(i)  General  densities  —  identified  through  fit  to  all  7  loops; 

(ii)  Nornral/lognormal  densities  —  identified  through  fit  to  all  7  loops; 

(iii)  General  densities  —  identified  through  fit  to  symmetric  major  loop  and  Rayleigh  loop; 

(iv)  Normal/lognormal  densities  —  identified  through  fit  to  major  loop  and  Rayleigh  loop. 
Convergence  was  achieved  in  all  cases  with  a  total  of  Nt  =  Nj  =  80  quadrature  points. 

General  Densities  —  Identification  Using  Full  Data  Set 

We  first  consider  the  model  performance  when  the  full  set  of  data  (all  7  loops)  was  used  to 
estimate  the  Ni  +  Nj  + 1  parameters  [rj,  u\ (EC1),  ■  ■  ■ ,  (ECn  ),  ^(-Eq),  •  •  • ,  V2(EiN. )]  through  solution 
of  the  nonlinear  constrained  optimization  problem  (24)  which  enforces  positivity  in  all  At,;  +  Nj  +  1 
parameters  and  symmetry  in  z^2.  The  resulting  model  fit  is  plotted  in  Figure  9  and  the  densities 
v\  and  v-2  are  shown  in  Figure  10.  It  is  observed  that  the  model  accurately  characterizes  both  the 
symmetric  major  loop  behavior  and  the  biased  minor  loop  responses  including  the  nested  biased 
minor  loop.  The  slight  oscillations  in  the  model  fit  preceding  saturation  result  from  the  optimization 
of  the  densities  to  accommodate  minor  loop  behavior  and  are  not  a  manifestation  of  numerical 
stability  for  large  stepsizes  in  the  input  fields.  For  later  comparison  with  fits  obtained  using  the 
lognormal/normal  densities,  the  residual  (39)  is  noted  to  be  77.  =  0.0057  as  summarized  in  Table  2. 

It  should  be  noted  that  when  traversing  from  the  endpoint  of  the  biased  minor  loop  to  the 
saturation  value,  no  a  priori  knowledge  of  the  connection  point  is  required  by  the  model.  The 
discontinuous  change  in  slope  is  automatically  incorporated  by  the  combination  of  energy-based 
hysterons  or  kernels  and  the  densities  used  to  incorporate  material  and  field  nonhomogeneities. 
The  automatic  enforcement  of  the  ‘deletion’  property  when  exiting  minor  loops  is  important  both 
for  material  characterization  and  the  development  of  model-based  compensators  for  linear  feedback 
control  design. 

Lognormal/Normal  Densities  —  Identification  Using  Full  Data  Set 

Using  the  full  set  of  data  from  the  7  loops,  we  also  estimated  the  parameters  77,  c,  b,  C  and 
Ec  in  the  discretized  model  (17)  in  which  zq  and  V2  are  specified  as  the  lognormal  density  (13) 


Electric  Field  (MV/m) 

Figure  9:  PZT5H  data  and  model  with  general  densities  v\  and  z/2  estimated  through  a  fit  to  full 
data  set. 
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(i) 

(ii) 

(iii) 

General  Densities  u\ ,  zz2  (A7; 

II 

II 

OO 

0.0057 

0.0115 

0.0050 

Lognormal/Normal  Densities  v\ ,  z^2 

0.0132 

0.0136 

0.0094 

General  Product  Density  u 

(. Ni  =  Nj  =  24) 

0.0028 

— 

— 

General  Product  Density  v 

(• Ni  =  Nj  =  48) 

0.0106 

— 

— 

General  Product  Density  v 

-  Regularized  (TVj  = 

-Nj- 

=  24) 

0.0037 

— 

— 

General  Product  Density  v 

-  Regularized  (TV*  = 

=  Nj  = 

=  48) 

0.0024 

— 

— 

Table  2:  Residuals  1Z  given  by  (39)  with  u\ ,  V2  estimated  using  the  nonlinear  constrained  minimization 
algorithm  (24)  and  v  identified  through  the  quadratic  programming  algorithm  (33)  and  regularized 
formulation  (38).  (i)  Identification  and  residuals  over  all  7  data  sets,  (ii)  Identification  using  major 
symmetric  and  Rayleigh  loops  —  residual  for  all  7  data  sets,  (iii)  Identification  using  major  symmetric 
and  Rayleigh  loops  —  residual  for  major  symmetric  and  Rayleigh  loops. 

and  normal  density  (15).  This  yields  the  parameter  values  r/  =  8.9  x  106,  Ec  =  7.6  x  105  V/m, 
c  =  0.237  V2/m2,  b  =  1.26  x  105  V2/m2,  C  =  1.4  x  10-12  and  resulting  model  fit  presented  in 
Figure  11.  The  corresponding  densities  are  plotted  in  Figure  12.  A  comparison  between  Figure  11 
and  Figure  9  illustrates  a  significant  loss  of  accuracy  when  characterizing  the  minor  loop  behavior 
with  this  a  priori  choice  of  densities.  This  is  quantified  by  the  residual  of  7 Z  =  0.0132  which  is  over 
twice  the  residual  obtained  with  general  densities. 

A  comparison  between  the  general  densities  plotted  in  Figure  10  and  the  lognormal/normal 
densities  plotted  in  Figure  12  indicates  the  same  tendency  observed  in  Section  4.1  for  PZT4  - 
the  two  sets  of  densities  exhibit  the  same  general  qualitative  behavior  but  the  general  densities  are 
significantly  less  regular  than  the  lognormal  and  normal  functions. 

General  Densities  —  Identification  Using  Symmetric  Major  and  Rayleigh  Loop  Data 

To  illustrate  the  predictive  properties  of  the  model  with  densities  estimated  from  limiting  data, 
we  consider  model  construction  using  data  from  the  symmetric  major  loop  and  Rayleigh  loop.  As 
illustrated  in  Figure  13,  the  resulting  model  fit  to  these  two  loops  is  highly  accurate  as  quantified  by 
the  residual  7 Z  =  0.0050  over  just  these  loops.  The  model  predictions  for  the  remaining  biased  minor 
loops  is  less  accurate,  however,  than  that  obtained  when  the  full  data  set  was  used  for  identification  — 
see  Figure  9.  This  loss  of  accuracy  is  quantified  by  the  total  residual  1Z  =  0.0115  which,  as  compiled 
in  Table  2,  is  over  twice  that  obtained  when  the  full  data  set  was  employed  for  identification. 


Figure  10:  Coercive  and  effective  field  densities  estimated  through  a  fit  to  full  data  set. 
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Figure  11:  PZT5H  data  and  model  prediction  with  lognormal  density  v\  and  normal  density  V2 
having  parameters  estimated  through  a  fit  to  full  data  set. 


Finally,  we  note  that  because  the  accuracy  of  biased  minor  loops  in  this  regime  is  highly  dependent 
on  the  accuracy  of  the  excursion  point,  even  minor  discrepancies  in  the  symmetric  loop  fit  can 
translate  into  large  minor  loop  errors  in  regions  having  a  large  gradient  ||I. 

Lognormal/ Normal  Densities  —  Identification  Using  Major  and  Rayleigh  Loop  Data 

In  the  final  case,  a  least  squares  fit  to  the  symmetric  major  and  Rayleigh  loop  data  was  used 
to  estimate  the  parameters  =  1.0  x  10',  Ec  =  7.6  x  105  V/m,  c  =  0.239  V2/m2,  b  =  1.22  x 
105  V2/m2,  C  =  1.4  x  10~12  in  the  discretized  model  (17)  constructed  using  the  lognormal  and 
normal  densities.  The  model  behavior,  plotted  in  Figure  14,  indicates  a  loss  of  accuracy  in  both  the 
fitted  loops  and  the  predictions  of  biased  minor  loop  behavior.  This  is  substantiated  by  the  residual 
values  of  0.0094  over  the  fitted  loops  and  0.0136  over  the  full  data  set.  Hence  while  the  model 
constructed  using  the  lognormal  and  normal  densities  is  sufficiently  accurate  for  many  applications, 
a  significant  improvement  in  accuracy  is  provided  by  the  identification  of  general  densities  using  more 
comprehensive  data  sets. 


Figure  12:  Lognormal  coercive  and  normal  effective  field  densities  estimated  through  a  fit  to  full 
data  set. 
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Electric  Field  (MV/m) 

Figure  13:  PZT5H  data  and  model  predictions  with  general  densities  v\  and  1/2  estimated  through  a 
fit  to  symmetric  major  and  Raleigh  loop  data. 

4.3  PZT5H  Characterization  —  Product  Density  v 

In  Section  3.3  we  summarized  a  solution  technique  for  estimating  the  product  density  v  based  on 
exploiting  the  linear  parameterization  to  reformulate  the  constrained  minimization  problem  as  the 
quadratic  programming  problem  (36).  This  offers  the  advantage  of  providing  a  well-established 
solution  framework  but  has  the  disadvantage  of  significantly  increasing  the  number  of  parameters  — 
Ni  ■  Nj  versus  Nj  +  Nj  +  1  —  and  can  be  prone  to  ill-conditioning. 

We  illustrate  in  this  example  the  use  of  this  formulation  to  estimate  v  in  (25)  using  all  7  loops 
of  the  PZT5H  data  originally  plotted  in  Figure  1.  To  assess  the  effects  of  conditioning,  we  consider 
two  discretization  levels  —  Np  =  Nq  =  6  (Nj  =  Nj  =  24)  and  Np  =  Nq  =  12  ( IV \  =  Nj  =  48)  - 
both  of  which  are  relatively  coarse.  The  model  fits  and  estimated  product  densities  for  the  two  cases 
are  plotted  in  Figure  15  and  16.  As  summarized  in  Table  2,  the  7-loop  residuals  for  Nj  =  Nj  =  24 


Figure  14:  PZT5H  data  and  model  prediction  with  lognormal  density  u\  and  normal  density  z/2 
having  parameters  estimated  through  a  fit  to  symmetric  major  and  Raleigh  loop  data. 
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(b) 


Figure  15:  PZT5H  data  and  model  prediction  with  general  product  density  v  estimated  using  the 
quadratic  programming  formulation  (34)  with  data  from  all  7  loops,  (a)  Nt  =  Nj  =  24,  and  (b) 
Ni  =  Nj  =  48. 


and  Ni  =  Nj  =  48  are  respectively  0.0028  and  0.0106.  The  fits  in  Figure  15  in  combination  with 
the  residual  values  lead  to  the  following  conclusions:  (i)  the  accuracy  with  the  coarse  discretiza¬ 
tion  Ni  =  Nj  =  24  is  comparable  to  that  obtained  through  solution  of  the  nonlinear  constrained 
optimization  problem  (24)  as  illustrated  in  Figure  9,  and  (ii)  the  accuracy  degrades  significantly  as 
discretization  levels  are  refined.  The  latter  reflects  the  inherent  ill-conditioning  associated  with  the 
algorithm. 

The  densities  plotted  in  Figure  16  demonstrate  the  validity  of  the  decay  assumption  (iii)  in  (20) 
and  restriction  of  the  solution  operator  to  compact  domains. 

Hence  while  this  approach  can  be  highly  accurate  and  relatively  efficient  for  certain  discretization 
levels,  care  must  be  taken  to  avoid  degradation  in  accuracy  which  can  occur  as  quadrature  levels 
and  associated  parameter  numbers  are  increased.  This  motivates  consideration  of  the  regularized 
formulation  (38)  which  is  illustrated  in  Section  4.4. 


Figure  16:  General  product  density  v  with  (a)  N%  =  Nj  =  24,  and  (b)  N%  =  Nj  =  48. 
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4.4  PZT5H  Characterization  —  Product  Density  v 

The  regularized  least  squares  formulation  (38)  stabilizes  the  pseudoinverse  associated  with  the  inverse 
problem  by  shifting  singular  values  away  from  the  origin.  The  model  fits  obtained  using  (38)  with 
a  =  5  x  1020  for  TV,  =  Nj  =  24  (TV  =  576)  and  TV,  =  Nj  =  48  (TV  =  2304),  using  data  from  all 
7  loops,  are  illustrated  in  Figure  17.  The  corresponding  residuals  7 Z  =  0.0037  for  ATt  =  Nj  =  24 
and  1Z  =  0.0024  for  TV,  =  Nj  =  48  demonstrate  improved  accuracy  as  TV  =  TV,  •  Nj  is  increased 
which  is  in  contrast  with  the  behavior  for  the  unregularized  algorithm.  Hence  it  is  observed  that 
regularization  in  this  manner  eliminates  the  instability  illustrated  in  Figure  15  for  increasing  TV  and 
yields  highly  accurate  model  representations.  Further  details  illustrating  the  formulation  and  effect 
of  regularization  techniques  for  integral  hysteresis  models  are  provided  in  [29]. 


(a)  (b) 

Figure  17:  PZT5H  data  and  model  prediction  with  general  product  density  v  estimated  using 
Tikhonov  regularization  with  data  from  all  7  loops,  (a)  TVj  =  Nj  =  24,  and  (b)  TVj  =  Nj  =  48. 


5  Concluding  Remarks 

In  this  paper,  we  have  presented  a  general  framework  for  quantifying  hysteresis  and  constitutive 
nonlinearities  inherent  to  ferroelectric  compounds  by  combining  aspects  of  the  homogenized  free 
energy  framework  developed  in  [32]  and  general  Preisach  models  employing  general  densities  or 
measures  [10,  24,  25,  28].  In  the  first  step  of  the  development,  Helmholtz  and  Gibbs  energy  relations 
are  constructed  at  the  lattice  level.  By  directly  enforcing  necessary  conditions  to  minimize  energy 
for  regimes  in  which  thermal  relaxation  mechanisms  are  negligible,  employing  Boltzmann  theory 
to  balance  the  relative  thermal  energy  and  Gibbs  energy,  we  then  construct  mesoscale  polarization 
relations.  These  relations  can  subsequently  be  extrapolated  to  provide  suitable  macroscopic  models 
for  homogeneous,  single  crystal  compounds.  For  nonhomogeneous,  polycrystalline  materials  with 
variable  effective  fields,  the  local  relations  provide  kernels  or  hysterons  which  serve  as  a  basis  with 
material  and  field  variations  incorporated  through  the  assumption  that  local  coercive  and  effective 
field  parameters  are  manifestations  of  underlying  stochastic  distributions  rather  than  constants.  In 
the  original  development  [32] ,  models  were  based  on  the  a  priori  assumption  of  lognormal  and  normal 
densities  for  the  coercive  and  effective  fields  whereas  we  employ  here  general  densities  which  are 
estimated  through  fits  to  experimental  data.  This  provides  the  models  with  the  flexibility  of  Preisach 
models  while  retaining  the  advantages  associated  with  the  energy  formulation  for  the  hysterons. 
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Form  a  theoretical  perspective,  the  models  incorporate  reversibility  and  incorporate  certain  rate 
and  temperature-dependencies  directly  into  the  polarization  kernel  or  hysteron.  Furthermore,  they 
guarantee  closure  of  biased  minor  loops  in  quasistatic  or  low  frequency  drive  regimes,  and  ensure  the 
‘deletion’  property,  but  permit  accommodation  and  after-effects.  As  detailed  in  [9]  in  the  context 
of  magnetic  materials,  several  of  these  properties  are  not  inherent  to  classical  Preisach  models  and 
require  extensions  that  diminish  the  efficiency  of  extended  Preisach  formulations  for  transducer  design 
and  model-based  control  design. 

The  extension  of  the  theory  in  [32]  through  formulation  in  terms  of  general  densities  rather  than 
a  priori  choices  of  lognormal  and  normal  densities  has  significant  ramifications  from  the  perspective 
of  both  material  characterization  and  control  design.  For  material  characterization,  formulation  in 
terms  of  general  densities  provides  additional  accuracy  for  certain  hard  materials  and  drive  regimes 
involving  multiple  biased  minor  loops  as  illustrated  in  the  context  of  PZT4  and  PZT5H  data.  For 
model-based  control  design,  the  advantages  are  even  more  profound  —  the  model  formulated  in  terms 
of  general  densities  exhibits  a  linear  dependence  on  the  parameters  whereas  the  model  constructed 
from  lognormal  and  normal  densities  exhibits  a  nonlinear  dependence  on  parameters.  In  the  latter 
case,  one  is  limited  to  nonlinear  adaptive  parameter  estimation  techniques  of  the  type  developed 
in  [14,  15,  22]  and  there  are  few  options  if  considering  adaptive  control  design.  For  the  linear 
parameterization  associated  with  the  general  density  formulation,  one  can  consider  a  broad  range  of 
adaptive  identification  and  control  designs  [13,  35]  which  extends  significantly  the  flexibility  of  the 
method. 

It  should  be  noted  that  whereas  the  identification  of  general  densities  using  constrained  optimiza¬ 
tion  algorithms  can  be  significantly  more  expensive  than  identification  of  parameters  in  the  lognormal 
and  normal  densities,  the  two  formulations  are  equally  efficient  to  implement  once  parameters  have 
been  identified.  Hence  both  formulations  are  viable  for  subsequent  use  in  material  characterization, 
transducer  and  system  design,  and  model-based  control  design  and  implementation. 
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